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ABSTRACT 


The  radar  cross  section  (RCS)  of  various  geometries  with 
surface  resistivity  is  synthesized  using  a  physical  optics 
approximation.  Using  the  RCS,  a  set  of  equations  are 
developed  for  reconstructing  the  resistivity  across  the 
surface.  The  reconstructed  resistivity  is  then  compared  with 
the  original  resistivity  and  the  RCS  recalculated  using  the 
synthesized  resistivity  to  validate  the  accuracy  of  the 
synthesis  procedure. 
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1.  INTRODUCTION 


During  the  middle  1950' s,  the  idea  of  making  platforms 
unobservable  or  less  observable  to  radar  was  initiated  by  Bill 
Bahret  of  Wright  Air  Development  Center.  [Ref.  1]  The  initial 
motivation  for  low-observability  came  from  the  need  to 
increase  screening  provided  by  a  jammer.  [Ref.  1] 
Subsequently,  a  great  deal  of  research  in  the  area  of  radar 
cross  section  (RCS)  began. 

There  are  four  common  methods  for  reducing  the  RCS  of  a 
given  platform.  They  are  [Ref.  2] : 

1.  shaping, 

2.  surface  material  selection, 

3.  active  cancellation,  and 

4.  passive  cancellation. 

Shaping  to  reduce  RCS  requires  reducing  the  surface  area  that 
is  perpendicular  to  the  radar  source.  Surface  material 
selection  requires  the  addition  of  lossy  materials  to  the 
surface  of  the  target.  Active  cancellation  is  a  method  of 
reproducing  the  signal  180*  out  of  phase  from  the  radar. 
Passive  cancellation  is  the  addition  of  a  secondary  scatterer 
to  cancel  the  reflected  field  from  the  primary  target. 

In  this  thesis,  a  synthesis  method  that  can  be  used  for 
surface  material  selection  is  addressed.  That  is,  given  the 
monostatic  RCS  of  a  platform,  a  set  of  equations  are  presented 
which  can  be  used  to  choose  the  required  material.  Selecting 
a  surface  material  that  reduces  the  RCS  of  a  platform  requires 
the  use  of  a  radar  absorbing  material  (RAM) .  Other  types  of 
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materials  used  in  the  reduction  of  RCS  include  resistive  and 
conductive  films,  described  in  terms  of  ohms  per  square  (Q/D) 
or  mhos  per  square  (U/D) ,  respectively. 

Metal  surfaces  such  as  those  used  for  aircraft  and  ships 
closely  approximate  a  perfect  electric  conductor  (PEC) .  The 
addition  of  RAM  to  a  surface  causes  it  to  have  an  inductive  or 
capacitive  surface  impedance  rather  than  that  of  a  PEC. 
Impedance  surfaces  can  absorb  radar  frequency  (RF)  energy 
thereby  causing  a  reduction  in  the  scattered  energy. 

An  impedance  surface  supports  both  equivalent  magnetic 
and  electric  currents.  However,  there  are  two  special  cases 
of  impedance  surfaces  that  only  support  either  equivalent 
magnetic  or  electric  currents:  a  resistive  surface  does  not 
support  equivalent  magnetic  currents,  and  a  conductive  sheet 
does  not  support  electric  currents.  The  resistive  surface  is 
the  focus  of  this  thesis,  because  its  use  greatly  simplifies 
the  equations  used  in  the  synthesis  procedure. 

A  synthesis  procedure  can  be  developed  based  on  the 
desired  RCS  for  a  particular  target.  The  surface  is  separated 
into  triangular  subsections,  or  patches.  The  surface 
impedance  across  each  subsection  is  assumed  to  be  constant  but 
can  vary  from  patch  to  patch.  A  matrix  equation  is  derived 
consisting  of  a  vector  of  the  known  scattered  field,  a  target 
surface  scattering  matrix,  and  a  surface  impedance  vector.  By 
solving  the  matrix  equation,  the  impedance  vector  is  obtained 
and  the  material  properties  of  the  surface  can  be 
determined. 

In  Chapter  II,  the  theoretical  background  required  to 
evaluate  the  RCS  of  a  surface  is  presented.  Surface 
impedances,  reflection  and  transmission  coefficients  and  the 
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physical  optics  approximation  are  examined. 

In  Chapter  III/  the  synthesis  procedure  is  developed. 
The  procedure  is  applied  to  doubly  curved  surfaces  such  as 
spheres  or  portions  of  spheres. 

In  Chapter  TV,  the  results  of  the  synthesis  procedure  are 
compared  to  the  exact  solution.  Finally,  in  Chapter  V,  an 
analysis  of  the  results  is  given  as  well  as  recommendations 
for  future  work  in  this  area. 
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n.  THEORETICAL  BACKGROUND 


In  this  chapter,  the  basic  equations  of  electromagnetic 
scattering  are  presented. 

A.  RCS  FROM  THE  SCATTERED  FIELD 

The  RCS  of  a  target  can  be  determined  once  its  scattered 
field  is  known  [Ref.  3] 

lim  (2.1) 

o  =  ( = - )  . 

|E.|^ 

The  scattered  field  arises  from  the  currents  induced  by  the 
plane  wave  on  the  target  surface.  The  scattered  electric 
field  far  from  the  target  is  given  in  terms  of  the  electric 
(7)  and  equivalent  magnetic  (M)  currents  by  [Ref.  4] 


ke 

Ee=-j^(L,+nN,)  , 

4jir 


(2.2) 


E*=j-^(Le-TlN^)  , 

4jir 


(2.3] 


where : 


N0= J f  (J^cos6cos(t>+JyCOs9sin4>-J^sm0)e-'^''*“'''ds'  , 


(2.4) 
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Lg= JJ"  (M^cosdcos4>+MyCOsdsm4>-M^sm6)e^  , 


(2.6) 


L^= / /eM^sm<t)+MyCos(t))e^'*“»ds'  .  (2.7) 

In  equations  (2.2)  through  (2.7)  k  is  the  wave  number  (2n/X, 
X  being  the  wavelength) .  The  geometry  is  shown  in  Figure  1 
below.  The  angles  0  and  ({)  are  the  standard  spherical  angles. 
For  the  remainder  of  this  thesis,  only  resistive  sheets  are 
considered  and  therefore,  X^=My=M2=0. 


Figure  1 .  Coordinate  system. 
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B.  SURFACE  IMPEDANCES 


Traditionally,  most  aircraft  and  ships  were  made  of 
materials  which  are  considered  to  be  perfect  conductors. 
Conducting  targets  generally  have  a  high  RCS  relative  to  those 
of  other  materials.  In  RCS  reduction  (RCSR) ,  the  selection  of 
materials  for  a  target  can  reduce  the  scattered  field.  One 
method  is  to  coat  a  perfectly  conducting  target  with  RAM. 
This  yields  a  heavy  target  and  is  not  an  effective  way  to 
reduce  RCS.  Another  approach  is  to  use  composite  materials 
with  some  loss  introduced.  The  latter  approach  is  more 
efficient,  since  composites  are  frequently  used  in  the 
manufacture  of  aircraft  and  ships.  Composites  can  be 
characterized  by  their  complex  index  of  refraction,  tv,  which 
is  related  to  the  permittivity  and  permeability  of  the 
material 


(2.8) 


where  e,  is  the  permittivity  of  the  material  and  (i,  is  the 
permeability  of  the  material.  The  permittivity  and 
permeability  are  usually  complex  niambers  and  it  is  the 
imaginary  parts  which  contribute  to  losses.  [Ref.  2]  Thin 
layers  of  lossy  material  can  be  applied  to  a  dielectric  to 
cause  the  surface  to  be  resistive  and,  no  magnetic  currents 
will  flow  on  the  surface. 

If  the  target  consists  of  composite  materials  then  there 
is  some  penetration  of  the  field  into  the  target.  The 
integrals  in  equations  (2.4)  and  (2.5)  become  volume  integrals 
because  the  currents  exist  below  the  surface.  This  is  very 
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inconvenient  computationally,  and  the  problem  can  be  avoided 
using  the  surface  impedance  approximation,  also  known  as  the 
Leontovich  boundary  condition 


(2.9) 


and  are  the  total  tangential  fields  on  the  outer 

surface  of  the  target  and  r\^  is  its  surface  impedance  in  ohms 
(ri.=0  if  the  surface  is  a  PEC).  [Ref.  3]  Thus  the  original 
target  surface  can  be  replaced  by  a  surface  with  impedance  Tj* 
and  the  same  scattered  fields  result.  Note  that  t),  depends  on 
the  location  on  the  surface  and  the  incidence  direction. 

Using  equation  (2.9),  a  set  of  vector  equations  can  be 
derived 

E-fi  (a  •  i)=Ti,7^  ,  (2.10) 

-  M.  (2.11) 

H-a  (a  •  H)^ — - . 

For  a  resistive  surface,  this  simplifies  to 

a  X  (a  X  e)=-rJ,  ,  (2.12) 

a  X  (a  X  H)=o  ,  (2.13) 

where  is  the  resistivity  of  the  material. 
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C.  REFLECTION  AND  TRANSMISSION  COEFFICIENTS 

The  reflection  coefficients  quantify  the  amount  of 
reflection  and  transmission  of  an  electromagnetic  wave  that 
occurs  at  an  interface  between  two  different  impedances.  For 
a  PEC,  the  incident  wave  is  totally  reflected,  yielding  a 
reflection  coefficient  of  -1.  The  reflection  (D , 

absorption  (A) ,  and  transmission  (T)  coefficients  must  obey  the 
conservation  of  energy  as  shown  below  [Ref.  2]: 

|r|^+|A|^+|T|^=l  .  (2.14) 

For  a  lossless  material,  A=0. 

The  transmission  and  reflection  coefficients  are  easily 
derived  with  the  assumption  that  the  surface  is  locally  flat. 
Two  principle  polarizations  must  be  considered:  1)  transverse 
magnetic  to  the  z  axis  (TM)  and  2)  transverse  electric  to  the 
z  axis  (TE)  .  An  arbitrary  polarized  wave  can  be  handled  as  a 
combination  of  the  two  principal  polarizations. 


1.  TM  Polarization 

For  TM  polarization,  the  incident  wave  is 


(2.15) 


where  k;  is  the  incident  wave  vector.  For  a  TM  polarized 
incident  wave,  the  reflection  coefficient  for  a  resistive  film 
is  [Ref.  5] 
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^TM 


-^ipCose. 

2R,+ti^cosej 


(2.16) 


where  0j  is  the  angle  of  incidence  with  respect  to  the  z  axis. 
If  no  absorption  occurs,  as  is  the  case  with  an  infinitely 
thin  sheet,  the  transmission  coefficient  is 


T  = 


2R. 


(2.17) 


2R, +11^0080. 


The  scattered  field  then  becomes 


(2.18) 


where  k  is  the  wave  vector  in  the  direction  of  observation, 


2.  TE  Polarization 

For  TE  polarization,  the  incident  wave  is 

(2.19) 

For  a  TE  polarized  incident  wave,  the  reflection  coefficient 
is  [Ref.  5] 


(2.20) 

2R  cos0.+T^ 

6  1  '  0 

With  no  absorption,  the  transmission  coefficient  is  given  by 
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(2.21) 


T 


TE 


2R  COS0. 

S  i 

2R,cos0j+ii^ 


and  the  scattered  field  becomes 


(2.22) 


D.  PHYSICAL  OPTICS  APPROXIMATION 

The  field  incident  on  an  illuminated  object  causes 
currents  to  flow  on  its  surface.  To  solve  exactly  for  the 
scattered  electric  field,  these  currents  must  be  )cnown.  It  is 
a  lengthy  and  complicated  process  to  solve  for  the  currents. 
However,  the  incident  field  and  the  boundary  conditions  are 
known.  Therefore  an  approximation  can  be  made  using  the  known 
incident  field.  The  standard  physical  optics  approximation  is 
given  by  [Ref.  3] : 


J,=2  fi  X  Hj  .  (2.23) 

Equation  (2.23)  holds  for  the  portion  of  the  object  which 
is  illuminated.  Physical  optics  approximation  assumes  the 
current  along  the  surface  drops  to  zero  on  areas  of  the  object 
which  are  shaded  from  the  illuminated  field.  In  actuality, 
currents  do  not  drop  abruptly  to  zero.  The  currents  slowly 
decay  towards  zero,  but  for  electrically  large  bodies  the 
error  is  negligible. 
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The  physical  optics  approximation  does  not  account  for 
waves  traveling  along  the  surface  of  the  object.  Therefore, 
there  is  a  decrease  in  RCS  as  the  incident  angle  moves  towards 
surface  grazing  and  the  sidelobes  of  the  RCS  pattern  are  lower 
than  the  actual  values.  The  specular  lobe  level  of  the  RCS 
pattern  for  flat  or  gently  curved  surfaces  is  generally  within 
a  few  percent  of  the  actual  value. 
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m.  FORMULATION  AND  SOLUTION 


The  analysis  of  scattering  from  an  arbitrary  finite 
doubly  curved  surface  requires  the  use  of  triangular 
subsections.  The  physical  optics  approximation  and  boundary 
conditions  are  used  to  obtain  the  synthesis  equations.  The 
synthesis  equations  are  then  verified  by  comparing  with 
calculated  fields.  The  geometries  considered  for  this 
calculation  are  a  flat  disk  and  a  spherical  cap.  The 
geometries  are  shown  in  Figures  2  and  3. 

The  incident  wave  is  assumed  to  be  a  planar  wave  with 
arbitrary  polarization.  Both  TM  and  TE  components  exist.  The 
incident  field  is  given  by 


E.=(0Ee‘+$E;)  e 


-jk,.r 


(3.1) 


where 

-jkj*r=jk(xsin6jCOS<J)j+ysin6jSin<J)j+zcos0j)  .  (3.2) 

Two  cases  are  considered  for  the  formulation  and 
solution:  bistatic  and  monostatic.  Bistatic  is  the  most 
general  case.  The  transmit  and  receive  antenna  for  the  radar 
are  not  located  at  the  same  position.  Monostatic  is  the  case 
in  which  the  transmit  and  receive  antenna  are  co-located. 
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For  each  case,  a  flat  disk  is  considered  first.  The  disk 
is  a  flat  surface  with  a  curved  perimeter  that  requires  the 
use  of  triangular  subsections  to  model.  Next,  the  more 
complex  shape  of  a  doubly  curved  surface  is  solved.  The 
solution  for  a  doubly  curved  surface  is  more  complex  because 
a  transformation  of  coordinates  is  required  for  each 
triangular  subsection.  This  transformation  gives  the  surface 
normal  in  the  transformed  z  axis  direction.  Also,  for  the 
doubly  curved  surface,  the  incident  wave  polarization  in  the 
individual  triangle  coordinate  systems  is  usually  not  pure  TM 
or  TE,  even  if  the  polarization  is  purely  TM  or  TE  in  the 
global  coordinate  system. 

A.  BISTATIC  CASE 
1.  Flat  Disk 

For  a  flat  disk  in  the  x-y  plane,  the  normal  is  always  in 
the  z  axis  direction.  Therefore,  the  boundary  condition  from 
equation  (2.9)  can  be  written  as: 

2:  X  (i  X  i)=-Rj,  (3.3) 
where  the  electric  field  in  the  above  equation  is  given  by 


Er(0T™E^$T^E;)e 


-jk,T 


(3.4) 


The  triple  products  in  equation  (3.3)  becomes 
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(3.5) 


t  TL  {1  X  E)=[(-jtcos0.sin<{)j+^cos0.cos<|)j)T^j^E0  + 
(-tcos0.sm<{>.-5^cos<|).cos0.)T^gE^‘]e  . 

By  substituting  Ttm  and  Tje  from  equations  (2.25)  and  (2.29) 
into  (3.5)  and  dividing  by  -R,,  equation  (3.3)  can  be  solved 
for  the  current  along  the  surface 


_  2(-icos0.cos4>i+^cos0.sm<J).)  = 

J  =  [ - i - ! - : — !-e;+ 

2R +n  COS0. 

K  *  O  1 

2(-tsin<|).cos0.-^cos0jSiii<().)  j 
- 

2R,cos0j+t)^ 


:3.6) 


The  two  principal  polarizations  are  considered  separately. 


0.  TM  Polarization 

For  TM  polarization,  the  <})  component  of  the  incident 
field  is  zero.  Therefore,  J,  in  equation  (3.6)  simplifies  to: 


_  2(-tcos0jCos<{>.+5^cos0jSiii(J).)  j  .^.7 

J  = - E^e  *  . 

2R  +11  COS0. 


For  use  in  equation  (2.4),  equation  (3.7)  is  separated  into  x 
and  y  components 


_  -2cos0.cos<t>.^  i  -j%.f 

2R,+fi^cos0. 

2cos0,sin({>.  J  _if., 
K - ^ - ^Ee‘e  ‘ 

^  2R^+t^i^cos0j 


(3.8) 


(3.9) 
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Equation  (2.4)  simplifies  to: 


rr  2cOS6.COS<t>jCOS0COS<()  : 

N  =  I  /  [ - J - i - Eg  + 

®  2R,+n„coste.  (3.10) 

2R +n  COS0. 

fi  '  O  1 

Portions  of  equation  (3.10)  are  not  dependant  on  surface 
location  and  can  be  removed  from  the  integral.  Equation 
(3.10)  then  reduces  to 


N0=2EgCOS0jCos0cos(<|).+<|>) 

(3.11) 

The  portion  of  equation  (3.11)  inside  of  the  integral  can  be 
denoted  as  I  for  simplicity 


r  p  jkr'cotilr 

f  - - ^ - dx'dy'  . 

•2R^+noCOS0j 


w  C 

2R^+Tl^cos0i 


dx^dy^  . 


(3.12) 


The  resistivity  of  the  material  is  restricted  to  be 
constant  on  each  triangle,  but  can  be  a  can  be  a  function  of 
location  on  the  surface  of  the  object.  Therefore,  it  must 
remain  within  the  integral.  Thus  the  resistivity  function  can 
be  defined  as  a  series  of  pulses  with  unknown  coefficients 
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(3.13) 


N 

0-1 

where  c  is  the  triangle  number  and  N  is  the  total  number  of 
triangles  on  the  surface  of  the  object.  The  coefficients 
are  the  unknown  expansion  coefficients  and  are  the  step 
functions 


|l.  x,,y„eS^ 
Ip,  otherwise 


(3.14) 


where  S,  is  the  area  of  triangle  c.  Substituting  equation 
(3.13)  into  (3.12)  gives 


^  jkr'cos^r 

c  c 

N 


dx^dy^  . 


ro^e-^nocose. 

0-1 


(3.15) 


Breaking  the  integral  into  a  sum  of  integrals  over  the 
individual  triangles  gives 


g  jkr  Wir 


dx'dy^s^I 

0-1 


0 


(3.16) 
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The  c*  term  of  the  sum  is 


llff _  _  (3.17) 

2r^+ti.cos0. 


The  denominator  of  equation  (3.17)  can  be  defined  as 


2r +T1  COS0. 
0*0  1 


(3.18) 


The  integrals  over  the  individual  triangles  can  be  reduced  to 
closed  form  as  described  by  Moreira  and  Prata.  [Ref.  6]  In 
general. 


I  =2S  - 5 - . 

D^(D,-Dp  D^(D,-Dp  D,D^ 

In  equation  (3.19),  the  C  and  D  coefficients  are  given  by  the 
following  equations: 


Co=A,.3 

(3.20) 

(3.21) 

(3.22) 

Drk(Bc.2-Bc.s) 

(3.23) 
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where  is  the  amplitude  value  of  the  integrand  of  (3.17)  at 
the  i*^^  vertex  of  triangle  c,  and  kBc^i  is  its  phase.  The 
quantities  are  defined  in  Figure  4.  The  primed  coordinates 
are  the  local  Cartesian  coordinates  with  z"  normal  to  the 
patch. 


Figure  4.  Integration  cell  geometry. 


Equation  (3.19)  has  singularities  which  can  be  removed  by 
using  Taylor  series  expansions  for  four  special  cases. 
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Case  1,  |Dp|<Lt  and  |Dq|>Lt: 


‘  'jD 


n»0 


n! 


n+l 


(3.24) 


Case  2,  (Dp|<Lt  and  |D,|<Lt: 


Ic=2S,e 
0  0 


'“•EE 

n»0  m-0 


aDp)“OD,r^ 

(n+m  +  l)!  " 


(3.25) 


Case  3,  |Dp|>Lt  and  ID,|<Lt: 


I.=2S,e''^-e^'E-^[-^G(n+l,-Dp]  . 
n-o  ii!  n+l 


(3.26) 


Case  4,  IDpl^Lt,  ID^ULt,  and  |Dp-Dq|<Lt: 


[j(D„-D„)]”  e^'C„ 

I.-2S.VE  -  ,  "  [-C.O(...D  . 

jD  --0  n!  (n+l) 


(3.271 


In  the  above  equations,  n  and  m  are  the  series  terms,  Lt  is 
a  small  number  that  determines  the  region  of  applicability  of 
the  series  approximation,  and  G  is  a  function  of  the  form: 


G(n,-w)=r^s“ej'*ds  . 

J  0 


(3.28) 
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For  the  calculations  in  this  thesis,  values  of  Lt-.005  and 
m=n=2  were  sufficient. 

Substitution  of  equations  (3.16)  and  (3.19)  into  (3.10) 
yields 

N  (3.29) 

Ng=2E9Cos0jCos0cos(4).-(t>)X,  aj,  . 

C-1 


From  equation  (2.2),  the  expression  for  the  co-polarized 
scattered  field  is 


-jkti  i 

^ - Ee*cos0.cos0cos(4).-<l))2r 

2tct  0-1 


(3.30) 


Finally,  using  I  in  place  of  the  sxam,  the  co-polarized  RCS  is 

(kil/  2 

Og0= - |cos0.cos0cos(4).-^)I|  .  (3.31) 

To  calculate  the  cross-polarized  scattered  field, 
equation  (2.5)  is  used  and  the  same  procedure  employed  to 
obtain  the  co-polarized  RCS  is  followed.  Equation  (2.5) 
reduces  to 


N^=2EgCos0.sin(<J)-4)j)I  . 

Equation  (2.3)  for  the  scattered  field  becomes 


(3.32) 


- E  0  CO  s0  jSm(<t>  -  <t)j)I 

2iir 


(3.33) 
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and  the  cross-polarized  RCS  is 


•  (3  34) 

b.  TE  Polarization 

For  TE  polarization,  the  0  component  of  the  incident 
field  is  zero.  Therefore,  J,  in  equation  (3.6)  simplifies  to 

_  2(-tsm(|)jCOS0.-tcose.sm<l).)^  i  (3.35) 

2R.c<,.e,.n.  ■ 

For  use  in  equation  (2.4),  (3.35)  is  separated  into  x  and  y 

components . 

-2(sm(l>jCose.)  i  .jt^.7 

J - -  9 

2R^cos0.+t|  ^ 

^  -2(co66.sm4>.)^  ■ 

^  2R,cose.-Hno 

Equation  (2.4)  simplifies  to 

r  r  ,  , 

N.=-2(cos0cos0.sin(4>-<{)i)  /  / - dx'dy'  . 

‘•^•'•2R,co60.+ti^ 


(3.36) 


(3.37) 
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The  equation  for  I  corresponding  to  equation  (3.14)  becomes 


g  jkr'oostg  -jV 

2R,cosej+Ti„ 


dx^dy^  . 


(3.39) 


Equation  (3.19)  for  Ij  remains  the  same.  However/  the 
denominator  is  not  the  same  and  therefore/  a  new  variable 
is  defined: 


1 


2R,cosej+ti, 


:3.40) 


The  equation  for  the  cross-polarized  scattered  field  becomes 


-jkr 


Ee;= 


2xr 


-E^‘cos0.cos0  sm(<|)  - 


(3.41) 


and  the  cross-polarized  RCS  is 


®  6<|> 


1  CO  s0  .CO  s6  sm(<|>  -  <t).)I  I  ^ 


(3.42) 


The  scattered  field  for  the  co-polarized  direction  is 
derived  in  a  similar  manner  yielding 


-jkn  j 

EA=-^-!-5 - E;cos0jcos(4>-«|)i)l 


'♦<!> 


Inr 


(3.43) 
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and 


(ktj  ,  (3.44) 

- |cos0.cos(<t>-<t>.)I|^  . 

% 

c.  Arbitrary  Polarization 

The  most  general  incident  polarization  is  a 
combination  of  TM  and  TE  components.  The  total  scattered 
field  in  the  0  direction  is 


e;=e3Ve 


f 


(3.45) 


or. 


-jkt)  e  j 

Eg  = - Eglcos0j[cos6cos(<j>-4)|)+sin(<|)-<{>j)]  . 

2jir 


(3.46) 


And  the  total  scattered  field  in  the  <j)  direction  is 


E  *=E  *  +E  ‘ 


(3.47) 


or. 


-jfcnpe 

1%T 


-jkr 

— E^‘lcos0j[cos0sin(4>-  4)j) +cos(4>  -  4‘i)]  • 


(3.48) 
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Equations  (3.47)  and  (3.48)  can  now  be  put  into  a  matrix  form. 


where 


[^00  ■^04>  ^0 


(3.49) 


-jkti  e'-’^  j 

A»g= - E0COs6;COsecos(<|)-<|)j)  , 

2jtr 


(3.50) 


-jkti  .e  i 

Ag.= - E^cos0.cos0sm(<{)-<{)j)  , 

2iir 


(3.51) 


-jkiie-^^  i 

-^*0= - Eecos0.sm(<t)-<t>j)  , 

2nr 


(3.52) 


-jkll  i 

- — - E^cos0.cos(4>-4>i)  , 


(3.53) 


s,-E  >,1. . 


(3.54) 
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and 


N 


S,-Eu  ■ 

C-1 


(3.55) 


The  resistivity  of  each  triangle  is  embedded  in  and  via 
equations  (3.54)  and  (3.55). 


2.  Doubly  Curved  Surfaces 

For  a  doubly  curved  surface,  the  surface  normal  vector  of 
each  triangular  subsection  can  be  different.  The  equations 
(3.19)  through  (3.27)  still  hold  if  the  local  triangle 
coordinates  are  used.  Therefore  a  transformation  of  the 
coordinates  must  be  performed  for  each  triangular  facet  on  the 
surface  from  the  global  system  (x,y,z)  to  the  local  system 
(x",y",z")  as  shown  in  Figure  5.  The  z"  axis  is  in  the  same 
direction  as  the  facet  normal  n^.  After  the  transformation  of 
coordinates,  the  scattering  calculations  are  made  for  each 
triangle,  the  fields  summed  over  all  triangles,  and  then  a 
transformation  is  made  back  to  the  original  (global) 
coordinates . 

The  incident  field  is  first  converted  to  global  Cartesian 
coordinates  [Ref.  4] 


Ex 

cos4>sm6 

COS0COS(t> 

-sin<j> 

e; 

siii6sin(j> 

cos0sm(t> 

COS(t> 

Ex 

COS0 

-sin0 

0 

(3.56) 
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z 


Figure  5.  Transformation  of  coordinates  from  a  global  (x,y,z)  to  a  local  (x",y",z") 
system  for  triangle  c. 


Next,  a  rotation  about  the  z  axis  is  performed  to  give  a  new 
coordinate  system  [x\y',z')  .  Then  a  rotation  about  the  y'  axis  is 
done  to  give  the  local-  coordinate  system  of  triangle  c  (x",y",z") 


x" 

cosa^  sina  0 

c  c 

cosp,  0  -sinp^ 

X 

y" 

= 

-sina,  cosa  0 

C  C 

0  1  0 

y 

o 

o 

_ 1 

sinp^  0  cosp^ 

z 
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The  angle  is  the  angle  of  rotation  about  the  z  axis,  and  the 
angle  p,  is  the  angle  of  rotation  about  the  y'  axis,  where 


a,=taix'^( 

C  ' 


_L_1) 

X  t 


(3.58) 


and 


p,=cos'‘(a^  X  i)  . 


(3.59) 


The  incident  field  is  also  transformed  to  the  local 
coordinate  system  using  the  same  matrices  as  in  equation 
(3.57).  After  all  transformations  are  complete,  equations 
(3.45)  through  (3.55)  are  applied.  Once  the  scattered  field 
for  each  triangle  is  calculated,  it  is  transformed  back  to  the 
global  coordinate  system  and  the  fields  from  all  triangles  are 
summed  vectorially.  Equations  (3.31),  (3.34),  (3.42),  and 

(3.45)  are  then  used  to  calculate  the  total  RCS  of  the 
surface . 

B.  MONOSTATIC  CASE 

For  a  monostatic  radar,  4)i=<}>  and  0i=0  and  the  scattered 
fields  reduce  to 


E  *= 
^0 


-jtn.e-*' 

2%T 


Eg‘lCOS^0 


(3.60) 
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and 


Inr 


-jkr 

— E^‘lcos0  . 


(3.61) 


In  equations  (3.60)  and  (3.61),  I  is  still  given  by  equations 
(3.15)  and  (3.38).  Recall  that  and  are  dependant  on  Oj 
and  location  on  the  surface.  In  order  to  synthesize  an 
equation  for  resistivity,  and  cannot  be  dependant  on  0;. 
This  implies  that  the  surface  composition  must  change  for  each 
incidence  angle  (ie.,  an  adaptive  target) .  One  approximation 
that  eliminates  the  dependence  on  0;  is  to  consider  cos0i=l. 
In  this  case. 


and  equations  (3.60)  and  (3.61)  reduce  to 


E  = 
^0  - 


2nr 


Ee*! 


(3.62) 


(3.63) 


2nr 


e;i 


(3.64) 


C.  SYNTHESIS  EQUATIONS 

In  applying  the  synthesis  procedure,  only  the  monostatic 
case  is  considered.  The  scattered  field  at  various  angles  is 
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known  from  equations  (3.63)  and  (3.64).  For  synthesis,  the 
scattered  fields  in  these  equations  are  known.  The  values  for 
I  and  subsequently  (=bj  are  to  be  determined.  Faros  has 
shown  that  in  providing  Eq  and  uniform  sampling  of  the 
field  in  direction  cosine  space  (DCS)  leads  to  less  numerical 
errors  than  sampling  in  9  and  cf).  [Ref.  5]  Therefore,  it  is 
more  convenient  to  work  in  DCS.  Figure  6  shows  the  domain  of 
DCS.  The  visible  region  is  the  area  within  the  unit  circle. 


v=sin0sm(p 


u=sin0cos(j) 


Figure  6.  Direction  cosine  space. 


To  generate  field  values  for  the  left  sides  of  equations 
(3.63)  and  (3.64),  DCS  is  divided  into  squares  of  area  du  by 
dv.  There  are  Nu  divisions  along  -the  u  axis  and  Nv  divisions 
along  the  v  axis.  The  center  of  each  square  is  the  sample 
point  for  calculating  the  scattered  field.  This  gives  a  total 
of  Nu  by  Nv  samples  of  the  scattered  field.  However,  only  the 
points  within  the  visible  region  are  valid  for  use  in  the 
synthesis  equations.  The  number  of  points  inside  the  visible 
region  is  defined  to  be  N.  The  angles  0  and  (j)  for  use  in 
calculating  the  scattered  field  are  found  from 

6  =  giii"^(u^+v^)  (3.65) 


4)=tan“^(— )  . 
u 


(3.66) 


The  un)cnowns  in  the  synthesis  problem  are  the  expansion 
coefficients,  a^.  There  is  one  value  of  for  each  triangular 
subsection.  Let  a  be  a  vector  of  N  elements  representing  N 
unlcnowns.  In  order  for  the  equation  to  be  solved  uniquely,  N 
unlcnowns  require  N  equations  and,  therefore,  the  number  of 
samples  of  the  scattered  field  must  also  be  N.  A  matrix 
ecjuation  follows  from  equations  (3.15),  (3.60),  and  (3.61): 
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(3.67) 


or 


where 


E=Xa 


E>E  •(e„,<{»,)=E  , 


(3.68) 


(3.69) 


=  [-^^(0Ea-‘l>E;)IJ|..„___^ 


(3.70) 


and  c,n=l,  2,  .  . .  ,N. 

The  X  matrix  is  square  and  can  be  inverted  and  used  to 
solve  for  the  a  matrix 


(3.71) 


Then  using  equation  (3.62),  each  element  in  the  a  matrix  is 
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used  to  solve  for  the  resistivity  over  each  triangle, 

If  1  1 

'’*"2  a  •'^2) 

0 

Equation  (3.72)  represents  the  synthesis  equation  for 
calculating  the  resistivity  of  an  object  given  the  required 
RCS. 
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IV.  COMPUTER  IMPLEMENTATION  AND  DATA  ANALYSIS 


This  chapter  discusses  the  computer  coding  of  the 
equations  derived  from  Chapters  II  and  III.  The  first 
sections  cover  the  calculation  of  the  RCS,  while  the  second 
half  addresses  the  synthesis  calculation. 

A.  RCS  FROM  THE  SCATTERED  FIELD 

Using  equation  (2.1),  the  RCS  can  be  calculated  from  a 
scattered  field.  The  incident  field  is  assirmed  TM  polarized 
and  the  monostatic  RCS  is  calculated  for  the  scattered  field 
obtained  using  equations  (3.60)  and  (3.61).  A  listing  of  the 
program  POTEST  is  included  in  the  appendix.  It  should  be 
noted  that  POTEST  calls  several  subroutines  which  are  also 
included  in  the  appendix. 

A  IX  X  IX  flat  PEC  plate  was  used  as  the  initial  test 
case.  The  patch  model  of  the  plate  is  shown  in  Figure  7.  In 
Figure  8,  the  approximate  (PO)  results  are  compared  to  the 
method  of  moments  (exact)  results.  The  specular  lobe  was 
within  a  few  dB,  but  the  sidelobes  are  slightly  larger  in  the 
exact  results.  This  is  due  to  the  fact  that  the  approximation 
used  in  equation  (3.3)  neglects  any  traveling  waves. 

The  second  test  was  to  place  some  resistivity  on  the 
surface  of  the  plate.  First,  a  constant  resistivity  of  377 
ohms  was  used.  A  comparison  of  the  PO  and  MOM  results  is  seen 
in  Figure  9.  Note  that  with  a  resistance  on  the  surface,  the 
sidelobe  differences  decrease. 
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Figure  7.  Flat  plate  in  the  x-y  plane. 


RCS  in  db 


The  final  test  was  to  increase  the  resistivity  linearly 
from  0  to  377  Q/U  from  the  center  of  the  plate  to  the  edges. 
This  was  done  using  the  equation 

R^=377(x>y.")  (4.1) 

where  andy^  are  the  midpoints  of  triangle  c.  Equation  (4.1) 
is  for  a  flat  surface  laying  in  the  x-y  plane,  centered  at  the 
origin.  The  results  are  shown  in  Figure  10. 

The  next  step  was  to  look  at  a  flat  disk  and  a  spherical 
cap  with  and  without  resistivity  to  establish  the  accuracy  of 
PO  relative  to  MOM.  The  models  used  for  a  flat  disk  and  a 
spherical  cap  are  shown  in  Figures  2  and  3,  respectively.  The 
results  for  a  flat  disk  are  shown  in  Figures  11  through  13. 
Figure  11  is  the  RCS  of  a  conducting  disk.  Figure  12  for  a 
disk  with  constant  resistivity  of  377  □/□,  and  Figure  13  is 
for  a  disk  with  resistivity  given  by  equation  (4.1)  .  Figures 
14  through  16  are  the  results  for  a  spherical  cap.  Figure  14 
shows  the  RCS  for  a  PEC  cap.  Figure  15  for  constant 
resistivity  of  377  Q/D  across  the  cap,  and  Figure  16  for 
linear  resistivity  using 


R  =377 
0 


_K_ 

^maz 


(4.2) 


where  is  the  arclength  from  the  pole  to  the  center  of 
triangle  c  and,  is  the  total  arclength  from  the  pole  to  the 
base  of  the  cap. 
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RCS  in  db 


RCS  in  db 


30 


Figure  11.  RCS  of  a  PEC  flat  disk  (TM  polarized  incident  wave,  monostatic). 
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RCS  in  db 


RCS  in  db 
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Figure  14.  RCS  of  a  PEC  spherical  cap  (TM  polarized  incident  wave,  monostatic). 
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RCS  in  db 


Figure  15.  RCS  for  a  spherical  cap  with  377  Q/n 
wave,  monostatic). 


resistivity  (TM  polarized  incident 


Figure  16.  RCS  of  a  spherical  cap  with  surface  resistivity  given  by  equation  (4.2)  (TM 
polarized  incident  wave,  monostatic). 
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B.  RCS  SYNTHESIS 

The  synthesis  computer  code  is  similar  to  the  program 
POTEST  with  changes  to  utilize  DCS  samples  rather  than  4)  and 
0.  The  increments  du  and  dv  were  calculated  so  that  N  would 
equal  the  number  of  triangles  on  the  surface  of  the  object. 
A  short  program,  NBOX,  was  used  to  insure  this.  After  the 
number  of  boxes  and  triangles  were  matched,  the  program  P0TC02 
was  used  to  calculate  the  field  scattered  from  the  object  at 
the  required  u  and  v  points.  The  scattered  field  was  then 
used  in  SYNTH  to  calculate  the  resistance  on  each  triangle 
based  on  the  given  field.  Listings  of  all  of  these  programs 
are  included  in  the  appendix. 

First,  a  flat  conducting  plate  was  again  used  as  a  test 

case.  For  this  calculation,  N=80,  du=dv=0.2.  The  synthesized 
resistivity  was  expected  to  be  0  Q/D  across  all  triangles,  and 
this  was  indeed  the  result.  Next,  the  plate  was  given  a 
constant  resistivity  of  377  Q/D.  The  synthesized  resistivity 
exactly  matched  the  expected  original  values.  Finally,  the 
plate  was  given  a  linear  resistivity  across  the  triangles. 
The  synthesized  resistivity  was  within  less  than  one  percent 
of  the  expected  result.  A  comparison  plot  of  the  original 
resistivity  and  the  synthesized  is  shown  in  Figure  17  on  the 
following  page.  Note  that  adjacent  triangles  have  different 
midpoints  and  therefore  the  resistivity  curve  oscillates. 
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Normalized  Resistivity 


Figure  17.  Original  resistivity  compared  with  synthesized  resistivity  (flat  plate,  surface 
resistivity  given  by  equation  (4.1)). 
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The  flat  disk  was  next  used  to  calculate  the  synthesized 
resistivity.  For  this  calculation,  N=156,  du=dv=0 . 1429.  With 
a  constant  resistivity  of  0  Q/D  and  377  fi/D,  the  synthesized 
resistivity  was  within  one  percent  of  the  original  value. 
This  is  shown  in  Figures  18  and  19  respectively.  With  a 
linear  resistivity  given  by  equation  (4.1),  the  synthesized 
resistivity  was  exactly  the  same  as  the  original  resistivity, 
as  shown  in  Figure  20. 

Finally,  resistivity  synthesis  was  conducted  for  a 
spherical  cap.  For  this  calculation,  N=156,  du=dv=0 . 1429 . 
The  synthesized  resistivities  for  cases  0  Q/D  and  377  Q/D 
were  within  one  percent  of  the  original  resistivities.  This 
is  shown  in  Figures  21  and  22.  The  results  for  a  linear 
resistivity  given  by  equation  (4.2)  were  exactly  the  same  as 
the  original  values.  This  is  shown  in  Figure  23. 

In  order  to  verify  the  synthesized  resistivity  would 
yield  the  correct  RCS,  a  comparison  was  made  between  the 
original  RCS  and  the  synthesized  RCS.  For  all  three 
geometries,  0  Q/D,  377  Q/D,  and  linear  resistivity,  the 
synthesized  resistivities  were  used  in  the  MOM  code  to  compute 
the  RCS  patterns.  The  results  are  shown  in  Figures  24  through 
32.  The  curves  labeled  Orig  are  the  PO  RCS  values  used  to 
initiate  the  synthesis.  Those  labeled  Synth  were  obtained 
using  the  synthesized  resistivity  in  the  MOM  code,  PATCH. 
These  figures  are  essentially  duplicates  of  Figures  8  through 
16  because  the  synthesized  resistivity  was  almost  identical  to 
the  starting  values. 
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Figure  18.  Original  resistivity  compared  with -synthesized  resistivity  (flat  disk,  0  Q/D 
resistivity). 
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Figure  19.  Original  resistivity  compared  with  synthesized  resistivity  (flat  disk,  377  Q/D 
resistivity). 
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Normalized  Resistivity 


Figure  20.  Original  resistivity  compared  with  synthesized  resistivity  (flat  disk,  surface 
resistivity  given  by  equation  (4. 1)). 
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Figure  21.  Original  resistivity  compared  with  S5mthesized  resistivity  (spherical  cap,  0  Q/D 
resistivity). 
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Normalized  Resistivity 


Figure  23.  Original  resistivity  compared  with  synthesized  resistivity  (spherical  cap,  surface 
resistivity  given  by  equation  (4.2)). 
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RCS  in  db 


RCS  in  db 


Figure  27.  Original  RCS  compared  with  synthesized  RCS  (flat  disk,  0  QfO  resistivity). 
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Figure  30.  Original  RCS  compared  with  synthesized  RCS  (spherical  cap,  0  Q/D  resistivity). 
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Figure  31.  Original  RCS  compared  with  synthesized  RCS  (spherical  cap,  377  Q/D 
resistivity). 
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RCS  in  db 


V.  CONCLUSIONS 


A  procedure  for  synthesizing  the  RCS  of  a  target  of  known 
shape  has  been  presented.  The  equations  were  derived  for  both 
flat  and  doubly  curved  surfaces  and  demonstrated  for 
monostatic,  and  purely  resistive  surfaces.  The  equations  were 
derived  using  the  physical  optics  approximation  and, 
therefore,  traveling  waves  were  neglected.  Furthermore, 
multiple  reflections  were  neglected. 

The  RCS  of  each  type  of  surface  was  first  calculated 
using  PO  and  then  compared  to  MOM.  Next,  the  RCS  of  each 
surface  was  calculated  in  DCS  so  the  scattered  field  could  be 
sampled  and  used  in  the  synthesis  equations.  The  synthesis 
equations  were  then  used  to  find  the  resistivity  across  the 
surface.  Finally,  the  resistivity  obtained  from  the  synthesis 
procedure  was  used  to  recalculate  the  RCS  and,  therefore, 
demonstrated  that  the  synthesis  procedure  was  accurate.  The 
synthesized  resistivities  were  in  excellent  agreement  with  the 
original  values  in  spite  of  all  of  the  approximations.  This 
is  because  the  chosen  resistivity  distributions  yielded  low 
RCS  sidelobes  and  did  not  support  traveling  waves.  Similarly, 
the  resistive  surface  transmits  more  of  the  incident  field 
than  it  reflects.  Therefore,  multiple  reflections  are  not 
important  and  can  usually  be  neglected  in  the  synthesis  of  low 
RCS  patterns. 

A  more  rigorous  solution  could  be  obtained  by  using  the 
MOM  as  the  basis  of  the  synthesis  equations.  [Ref.  3]  This 
would  require  solving  exactly  for  the  current  along  the 
surface  as  well  as  for  the  resistivity.  Future  work  in  this 
area  should  include  a  synthesis  procedure  based  on  MOM. 
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APPENDIX.  COMPUTER  CODES 


This  appendix  contains  all  of  the  computer  codes  used  for 
calculations  in  this  thesis.  The  code  POTEST. M  calculates 
the  RCS  of  a  surface  with  specified  resistivity.  The  code 
P0TC02.M  calculates  the  RCS  of  a  surface  in  DCS.  The  code 
BUILD. M  generates  the  geometry  of  the  surface.  The  code 
FACT.M  is  a  subroutine  used  to  calculate  the  factorial  of  a 
number  and  is  used  with  POTEST. M  and  P0TC02.M.  The  code  GF.M 
is  a  subroutine  used  to  calculate  equation  (3.28)  and  is  used 
with  POTEST. M  and  P0TC02.M.  The  code  RES.M  is  a  subroutine 
used  in  the  calculation  of  linear  resistivity  across  a  surface 
using  equations  (4.1)  and  (4.2)  and  is  used  with  POTEST. M  and 
P0TC02.M.  The  code  SYNTH. M  reconstructs  the  resistivity  using 
the  scattered  field  from  P0TC02.M.  The  code  NBOX.M  calculates 
the  number  of  partitions  of  DCS  that  are  within  the  visible 
region  for  a  given  du(=dv).  This  was  used  to  verify  that  the 
number  of  known  scattered  field  points  would  be  equal  to  the 
number  of  triangular  subsections  on  the  surface.  The  code 
POS.M  compares  the  original  RCS  with  the  synthesized  RCS. 
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%  POTEST. M 

%  PO  scattering  from  triangles 
clear  all 

%  illiamination  flag:  =0  external  face  only 
iflag=l; 

Lt=.001; 

Nt=2  ; 
wave=l; 
bk=2*pi/wave; 
rad=pi/180/ 

%  Incindent  wave  polarization 
Et=l+j*0;  %TM-z 

Ep=0+j*0;  %TE-z 

phi=0; 

phr=phi*rad/ 

cp=cos (phr) ; sp=sin (phr) ; 

%  Wave  ampltidue  at  all  vertices 

Co=l . ; 

build 

ntria=nfaces; 

nvert=nverts; 

x=xpts; 

y=ypts; 

z=zpts; 

title (' triangular  model  of  the  scattering  surface') 

xlabel ( 'x' ) 

ylabel ( *  y ' ) 

zlabel('z') 

for  i=l:nvert 

text (x (i) -max (x) /20, y (i) -max (x) 720, z (i) , num2str (i) ) 
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end 

hold  off 

disp('hit  return  to  continue  calculation') 
pause 

%  Define  position  vectors  to  vertices 
for  i=l:nvert; 

r  (i,  : )  =  [x  (i)  y  (i)  z  (i)  ]  ; 

end 

%Use  res.iti  program  to  calculate  resistivity 
res 

%  Get  resistivity,  edge  vectors  and  normals  from  edge  cross 
products 

for  i=l:ntria 
rs  (i)=0; 

%  Resistivity  function  for  a  flat  surface 
%rs  {i)=rss*sqrt  (xm(i)  ''2+ym(i)  ^2)7377/ 

%  Resistivity  function  for  a  doubly  curved  surface 
%rs  (i)=rss*arc  (i)V377; 

A=r (vind (i, 2) , : ) -r (vind (i, 1)  ,  : ) ; 

B=r (vind(i, 3) , : ) -r (vind (i, 1) , : ) ; 

C=r (vind (i, 3) , : ) -r (vind (i, 2) ,  : )  ; 

N (i,  : ) =cross (A, B) ; 

%  check  for  proper  direction  of  normal  vector 
if  N(i,3)<0 

N(i,  :)=-N(i,  :)  ; 

end 

%  Edge  lengths  for  triangle  "i" 
d(i, 1) =norm(A) ; 
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d (i, 2) =norm (B) ; 
d(i, 3) =norm{C) ; 
ss=.5*sum(d(i, : )  )  ; 

Area (i) =sqrt (ss* (ss-d (i, 1) ) * (ss-d (i,  2) ) * (ss-d (i,  3)  )  )  ; 
Nn=norm(N (i, : ) ) ; 

N(i,  :  )=N(i,  : )  /Nn; 

%transform  triangles  into  local  coordinates 
xl=[l  00]; 
yl=[0  1  0]; 
zl=[0  0  1]; 

end  %  end  of  triangle  loop 

%  Pattern  Loop 
start=0; 
stop=180; 
del=l; 

it=floor ( (stop-start) /del) +1; 
for  i=l:it 

theta (i) =start+ (i-1) *del; 
thr=theta (i) *rad; 
st=sin(thr);  ct=cos(thr); 
u=st*cp;  v=st*sp;  w=ct/ 
up=ct*cp;  vp=ct*sp;  wp=-st; 

%  Spherical  coordinate  system  unit  vectors 
R(i,  :  )  =  [u  V  w]  ; 

%  Change  Efield  into  global  rectangular  coordinates 
strtra (1, : ) = [st*cp  ct*cp  -sp]; 
strtra (2, : ) = [st*sp  ct*sp  cp]; 
strtra (3, :)  =  [ct  -st  0]  ; 
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sei= [0  Et  Ep] ; 
ei=strtra*sei ' ; 

%  correct  orientation  for  normal  after  90deg 
if  theta  (i)  ===91 
N=-N; 

end 

%  Begin  loop  over  triangles 
siimt=0; 
sump=0; 
for  m=l:ntria 
t (m, : ) =cross ( zl , N (m, : ) ) ; 
beta (ra) =acos (N (m, 3 ) ) ; 

alpha (m) =atan2 (t (m, 2) , t (m, 1) +le-6) -pi/2; 

%  Test  to  see  if  front  face  is  illuminated 
ndotk=N{m, : ) *R(i, : ) ' ; 
if  ndotk  >=  0  &  iflag==0  1  iflag~=0 
%translation  into  local  coordinates 
%first  rotate  about  z  axis 

tral (1, :)= [cos (alpha (m) )  sin (alpha (m) )  0]; 

tral (2/ : ) = [-sin (alpha (m) )  cos (alpha (ra) )  0] ; 
tral (3,:)=[0  0  1]; 

%  Rotation  about  yp  axis 

tra2 (1, : )=[cos (beta(m) )  0  -sin (beta (m) )] ; 
tra2  (2,  :  )  =  [0  10]; 

tra2(3,:)=[  sin (beta (m))  0  cos (beta (m) )] ; 
Dl=tral*R (i, : ) ' ; 

D2=tra2*Dl; 
upp=D2 ( 1 ) ; 
vpp=D2 ( 2 ) ; 
wpp=D2 ( 3 ) ; 


73 


%  Find  spherical  angles  in  local  coordinates 
stpp=sqrt  (upp''2+vpp''2)  *sign  (wpp)  ; 
ctpp=sqrt  (l-stpp''2)  ; 
thpp=acos (ctpp) ; 
phipp=atan2 (vpp, upp+le-10) ; 

%  Phase  at  the  vertices  of  triangle  m;  monostatic  RCS  needs 
"2" 

Dp=2*bk* ( (x (vind (m, 1 ) ) -x (vind (m, 3) ) ) *u+ . . . 

(y (vind (m, 1 ) ) -y (vind (m, 3) ) ) *v+ . . . 

(z (vind(m, 1) ) -z (vind (m, 3) ) ) *w) / 

Dq=2*bk* ( (x (vind (ra, 2) ) -x (vind (m, 3) ) ) *u+ . . . 

(y (vind(m/ 2) ) -y (vind(m, 3) ) ) *v+. . . 

(z (vind (m, 2) ) -z (vind (m, 3) ) ) *w) ; 

Do=2*bk* (x (vind (m, 3) ) *u+y (vind (m, 3) ) *v+z (vind (m, 3) ) *w)  ; 
%incident  field  in  local  rectangular  coordinates 
eip=tral*ei; 
eipp=tra2*eip; 

%  translation  of  E-field  into  local  spherical  coordinates 
strtra2= [stpp*cos (phipp)  ctpp*cos (phipp)  -sin(phipp); 

stpp*sin (phipp)  ctpp*sin (phipp)  cos (phipp) ;  ctpp  -stpp  0]; 
epp=strtra2^ (-1) *eipp; 

%Reflection  coefficients 
gamp=0; 

if  (2*rs (m) +ctpp) ~=0 
gamp=l/ (2*rs (m) +ctpp) ; 
end 
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gamn=l/ (2*rs (m) *cos (thpp) +1) ; 

%Surface  Currents 

Jxpp  (i,m)  =  (epp  (2)  *gamp-epp  (3)  *gamn)  *ctpp''2; 

Jypp  (i,m)  =  (epp  (2)  *gainp+epp  (3)  *garan)  *ctpp; 

%  Area  integral  for  general  case 

DD=Dq-Dp ; 

expDo=exp ( j  *  Do )  ; 

expDp=exp ( j  *Dp )  ; 

expDq=exp ( j  *Dq)  ; 

%  Special  Case  1 
if  abs (Dp)  <  Lt  &  abs (Dq) >=  Lt 
sic=0 . ; 
for  n=0:Nt 

sic=sic+  ( j*Dp)  ''n/fact  (n)  *  (-Co/  (n+l)  +expDq*  (Co*gf  (n,  -Dq)  )  )  ; 
end 

Ic  (i/iti)  =sic*2*Area  (m)  *expDo/j/Dq/ 

%  Special  Case  2 

elseif  abs (Dp)  <  Lt  &  abs(Dq)  <  Lt 
sic=0 . ; 
for  n=0:Nt 

for  nn=0:Nt 

sic=sic+ ( j*Dp) ^n* ( j*Dq) ^nn/fact (nn+n+2) *Co; 

end 

end 

Ic  (i/in)  =sic*2*Area  (m)  *expDo; 

%  Special  Case  3 

elseif  abs (Dp)  >=  Lt  &  abs(Dq)  <  Lt 
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sic=0 . ; 
for  n=0:Nt 

sic=sic+  ( j*Dq)  ''n/fact  (n)  *Co*gf  (n+1,  -Dp)  /  (n+l)  ; 
end 

Ic (i,m) =sic*2*Area (m) *expDo*expDp; 

%  Special  Case  4 

elseif  abs(Dp)  >=  Lt  &  abs(Dq)  >=  Lt  &  abs (DD)  <  Lt 
sic=0. ; 
for  n=0:Nt 

sic  =  sic+ ( j  *DD)  ^'n/fact  (n)  *  (-Co*gf  (n,  Dq) +expDq*Co/ (n+l )  )  ; 
end 

Ic (i,m) =sic*2*Area (m) *expDo/ j /Dq; 
else 

Ic  (i,m)  =2*Area  (m)  *expDo*  {expDp*Co/Dp/DD-expDq*Co/Dq/DD-Co/Dp 
/Dq)  ; 

end  %  end  of  special  cases 

test 

%Scattered  field  in  local  rectangular  coordinates 
Esxpp=Jxpp  (i,in)  *Ic  (i,m)  ; 

Esypp=Jypp (i,m) *Ic (i,m) ; 

Eszpp=0; 

%retranslation  to  global  rectangular  coordinates 
Espp (m, : ) = [Esxpp  Esypp  Eszpp] ; 

Esl=tra2''  (-1)  *Espp  (m,  : )  '  ; 

Es=tral''  (-1)  *Esl; 

%retranslation  to  global  spherical  coordinates 
%seipp  (m,  : )  =strtra''  (-1 )  *Es; 
seipp (m, 1) =0; 
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seipp (m, 2) =ct*cp*Es (1) +ct*sp*Es (2) -st*Es  (3) ; 
seipp (m, 3) =-sp*Es ( 1 ) +cp*Es (2 ) ; 

%siinv  over  all  triangles  to  get  scattered  field 
sumt=sumt+seipp (m, 2) ; 
sump=sump+seipp (m, 3) ; 

end  %end  of  illumination  test 
end  lend  of  triangle  loop 
Eth=sumt ; 

Sth(i)=10*logl0 (4*pi*abs (Eth) ^2+l.e-lO) ; 

Eph=sump; 

Sph (i) =10*logl0 (4*pi*abs (Eph) ^2+1 .e-10) ; 
end 

Sth  ( : )  =max  (Sth  ( : )  , -50)  ; 

Sph  ( : )  =max  (Sph  ( : )  ,  -50)  ; 

%plot  res 
figure (2) 

plot (theta, Sth, theta, Sph, ' — ' ) , grid 

title (' solid:  co-polarization  dashed:  cross-polarization') 
xlabel ( 'monostatic  angle,  theta  (deg)') 
ylabel  ( '  sigma/wavl''2  (dB)  ' ) 
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%  P0TC02.M 

%  PO  scattering  from  triangles 
clear  all 
format  long 

%  illumination  flag:  =0  external  face  only 
iflag=0; 

Lt=.001; 

Nt=2  ; 
wave=l ; 
bk=2*pi/wave; 
rad=pi/180; 

%  Incindent  wave  polarization 

Et=l+j*0;  %TM-z 

Ep=0+j*0;  %TE-z 

%  Wave  ampltidue  at  all  vertices 

Co=l . ; 

build 

ntria=nfaces; 
nvert=nvert s ; 
x=xpts; 
y=ypts; 
z=zpts; 

title (' triangular  model  of  the  scattering  surface') 

xlabel { 'x' ) 

ylabel ( 'y' ) 

zlabel  ( ' z' ) 

for  i=l:nvert 

text (x (i) -max (x) /20, y (i) -max (x) /20, z (i) ,num2str (i) ) 

end 

hold  off 
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dispChit  return  to  continue  calculation') 
pause 

%  Define  position  vectors  to  vertices 
for  i=l;nvert; 

r(i,:)  =  [x(i)  y(i)  z(i)]; 

end 

%Use  res.m  program  to  calculate  resistivity 
%res 

%  Get  resistivity/  edge  vectors  and  normals  from  edge  cross 

products 
for  i=l:ntria 
rs  (i) =0; 

%resistivity  for  a  plate 

%rs  (i)=rss*sqrt  ((xm(i)-.5)''2+(ym{i)-.5)''2)/377; 

%  resistivity  for  a  disk 

%rs  (i)=rss*sqrt  (xm (i)  "'2+ym (i)  ^^2 )  /  377 ; 

%  resistivity  for  a  cap 
%rs (i)=rss*arc (i) /377; 
ac (i)=l/ (2*rs (i) +1) ; 

A=r (vind(i, 2) , : ) -r (vind(i,  1)  ,  : )  ; 

B=r (vind(i, 3) , : ) -r {vind(i,  1)  ,  : )  ; 

C=r (vind (i/ 3) , : ) -r {vind(i, 2)  ,  : ) ; 

N(i, : )=cross (A,B) ; 

%  check  for  proper  direction  of  normal  vector 
if  N(i,3)<0 

N(i,  : )  =-N{i,  : )  ; 

end 

%  Edge  lengths  for  triangle  "i" 
d(i, l)=norm(A) ; 
d(i/  2)  =norm(B)  ; 
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d(i,  3)=norm(C) ; 
ss=.5*sum(cl(i,  : )  )  ; 

Area (i ) =sqrt (ss* (ss-d(i/ 1) ) * (ss-d(i,2) ) * (ss-d(i,3) ) ) ; 
Nn=noriii(N  (i,  : )  )  ; 

N(i,  :  )=N(i,  : )  /Nn; 

%transform  triangles  into  local  coordinates 

xl=[l  0  0]  ; 

yl=[0  1  0]  ; 

zl=[0  0  1]  ; 

end 

%  Pattern  Loop 
it=0; 

load  du.m; 
k=0; 

for  ul=-l : du: 1-du 

for  vl=-l :du: 1-du; 

it=it+l; 

Eth(it)=0; 

Eph(it)=0; 

%test  to  see  if  center  of  box  is  within  visible  region 
test=  (  (ul+du/2)  ''2+  (vl+du/2)  ''2)  ; 
if  test<=l 
k=k+l; 

phi=atan2 ( (vl+du/2) ,  (ul+du/2) )  ; 
cp=cos (phi) ;  sp=sin (phi) ; 
theta=asin (sqrt (test) ) ; 
st=sin (theta) ;  ct=cos (theta) ; 

%  Center  of  box  is  used  so  u  and  v  are  redefined  using  center 
of  %  box 

u=st*cp;  v=st*sp;  w=ct; 
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up=ct*cp;  vp=ct*sp;  wp=-st; 

%  Spherical  coordinate  system  unit  vectors 
R(it,  :  )  =  [u  V  w]  ; 

%  Change  Efield  into  global  rectangular  coordinates 

strtrad,  :  )  =  [st*cp  ct*cp  -sp]; 

strtra(2, : )=tst*sp  ct*sp  cp]; 

strtraO,  :  )  =  [ct  -st  0]; 

sei=  [0  Et  Ep] ; 

ei=strtra*sei ' ; 

%  Begin  loop  over  triangles 
for  m=l:ntria 
t  (m,  : )=cross (zl,N (m, : ) )  ; 
beta (m) =acos (N (m, 3 ) ) ; 

alpha (m) =atan2 (t (m,  2) , t (m, 1) +le-6) -pi/2; 
lc{k,m)=0; 

%  Test  to  see  if  front  face  is  illuminated 
ndotk=N(m, : ) *R(it, : ) ' ; 
if  ndotk  >=  0  &  iflag==0  1  iflag~=0 
%translation  into  local  coordinates 
% first  rotate  about  z  axis 

tral (1, : )  =  [cos (alpha (m) )  sin  (alpha (m) )  0]; 

tral (2,  :)  =  [-sin(alpha(m) )  cos (alpha (m) )  0]  ; 

tral (3, :)=[0  0  1] ; 

%  Rotation  about  yp  axis 

tra2 (1, :)= [cos (beta (m) )  0  -sin (beta (m) )] ; 

tra2 (2, : ) = [0  1  0] ; 

tra2(3,:)  =  [  sin  (beta  (m))  0  cos  (beta  (m)  )] ; 

Dl=tral*R (it-,  : )  ' ; 

D2=tra2*Dl; 
upp=D2 ( 1 ) ; 
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vpp=D2 { 2 ) ; 
wpp=D2.(3)  ; 

%  Find  spherical  angles  in  local  coordinates 
stpp=sqrt  (upp''2+vpp"2)  *sign(wpp)  ; 
ctpp=sqrt (l-stpp^2)  ; 
thpp=acos (ctpp) ; 
phipp=atan2 (vpp,upp+le-10)  ; 

%  Phase  at  the  three  vertices  of  triangle  m;  monostatic  RCS 
needs  %  "2" 

Dp=2*bk* (  (X (vind (m, 1 ) ) -x { vind (m, 3 ) ) ) *u+ . . . 

(y (vind (m, 1 ) ) -y (vind (m, 3) ) ) *v+ , . . 

(z(vind(m,l) ) -z (vind (m, 3) ) ) *w) ; 

Dq=2*bk* ( (x (vind (m, 2 ) ) -x (vind (m, 3 ) ) ) *u+ . . . 

( y (vind (m, 2 ) ) -y (vind (m, 3 ) ) ) * v+ . . . 

(z (vind(m, 2) ) -z (vind (m,  3)  )  )  *w)  ; 

1)0=2 *bk*  (x  (vind  (m,  3 )  )  *u+y  (vind  (m/  3 )  )  *v+z  (vind  (m,  3 )  )  w)  , 


%incident  field  in  local  rectangular  coordinates 

eip=tral*ei; 

eipp=tra2*eip; 

%  translation  of  E-field  into  local  spherical  coordinates 
strtra2= [stpp*cos (phipp)  ctpp*cos (phipp)  -sin(phipp)/ 

stpp*sin(phipp)  ctpp* sin (phipp)  cos (phipp);  ctpp  -stpp  0] ; 
epp=strtra2" (-1) *eipp; 

%  Area  integral  for  general  case 
DD=Dq-Dp; 
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expDo=exp ( j *Do) ; 
expDp=exp ( j  *Dp) ; 
expDq=exp ( j  *Dq)  ; 

%  Special  Case  1 
if  abs(Dp)  <  Lt  &  abs(Dq)>=  Lt 
sic=0 . ; 
for  n=0:Nt 

sic=sic+ ( j*Dp) ^n/fact (n) * (-Co/ (n+1) +expDq* (Co*gf  (n,  -Dq) ) )  ; 
end 

Ic (k,m) =sic*2*Area (m) *expDo/j/Dq; 

%  Special  Case  2 

elseif  abs(Dp)  <  Lt  &  abs(Dq)  <  Lt 
sic=0 . ; 
for  n=0:Nt 

for  nn=0:Nt 

sic=sic+  {j*Dp)  ''n*  (j*Dq)  ''nn/fact  (nn+n+2)  *Co; 

end 

end 

Ic (k,m)=sic*2*Area (m) *expDo; 

%  Special  Case  3 

elseif  abs(Dp)  >=  Lt  &  abs(Dq)  <  Lt 
sic=0 . ; 
for  n=0:Nt 

sic=sic+  (j*Dq)  ^n/fact  (n)  *Co*gf  (n+.l,-Dp)  /  (n+1)  ; 
end 

Ic (k, m) =sic*2*Area (m) *expDo*expDp; 

%  Special  Case  4 
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elseif  abs (Dp)  >=  Lt  &  abs (Dq)  >—  Lt  &  abs (DD)  <  Lt 
sic=0 . ; 
for  n=0:Nt 

sic  = 

sic+  ( j  *DD)  ''n/fact  (n)  *  (-Co*gf  (n,  Dq)  +expDq*Co/  (n+1 )  )  ; 
end 

Ic  (k,ni)  =sic*2*Area  (m)  *expDo/j/Dq; 
else 

Ic (k,m) =2*Area (m) *expDo* (expDp*Co/Dp/DD-expDq*Co/Dq/DD-Co/Dp 
/Dq)  ; 

end 

%  end  of  special  cases  test 
end  %end  of  illximination  test 
end  %end  of  triangle  loop 

%sum  over  all  triangles  to  get  scattered  field 
Eth(it)=Ic (k, : ) *ac'*Et; 

Eph (it) =Ic (k, : ) *ac' *Ep; 

%  correction  for  rounding  errors 
%Eth2 (it) =round (Eth (it) *le4) ; 

%Eth(it)=Eth2 (it)*le-4; 

end  %  end  of  test  for  visible  region 

Sth(it)=10*logl0  (4*pi*abs  (Eth (it)  /  (Et+le-10)  )  ''2+l.e-lO)  ; 

Sph (it)  =10* logic  (4*pi*abs  (Eph (it)  /  (Ep+l.e-10) )  ''2+l.e-lO)  ; 
end  %  end  of  v  loop 
end  %  end  of  u  loop 
Sth(:)=max(Sth(:),-50)  ; 

Sph  ( : )  =inax  (Sph  ( : )  ,  -50)  ; 

Sth=Sth'; 

Sph=Sph ' ; 

%plot  res 
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u=-l : du: 1-du; 
v=-l :du: 1-du; 

T=ones (1,  2/ du) ; 
nrow=2/du; 
ncol=2/du; 
for  k=l:ncol; 

kl=nrow* ( k- 1 ) + 1 ; 
k2=kl+nrow-l; 

X2(:,k)  =  (u(k)+du/2)*T'; 
Y2{:,k)  =  (du/2.+v)  '  ; 

Z2  ( : ,  k)  =Sth  (kl :  k2)  ; 

end 

figure (2) 
orient  tall 
surf (X2, Y2, Z2) 
view (-25, 55) 
colormap (hot) 
xlabel ( 'u' ) 
ylabel ( ' v' ) 
zlabeK'RCS  in  Db') 

Eth=Eth' ; 

Eph=Eph ' ; 
rs=rs ' ; 

Ethr=real (Eth) ; 

Ethi=imag (Eth) ; 

Ephr=real (Eph) ; 

Ephi=imag (Eph) ; 

Ici=imag (Ic)  ; 

save  Ethr.m  Ethr  -ascii  -double 
save  Ethi.m  Ethi  -ascii  -double 
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save  Ephr.m  Ephr  -ascii 
save  Ephi.m  Ephi  -ascii 
save  rs.m  rs  -ascii 
save  Icr.m  Ic  -ascii 
save  Ici.m  Ici  -ascii 
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%  BUILD. M 


%  Program  to  read  transposed  "buildnS"  data  and  plot  it  in 
MATLAB 

%  node  (x,y, z)  coordinates 
%  legend: 
load  xpts.m 
load  ypts.m 
load  zpts.m 
nverts=length (xpts) ; 

%  node  connection  list 
load  endl.m 
load  end2.m 
nedges=length (endl) ; 

%  face  connection  list 
load  nodel .m 
load  node2 .m 
load  node3.m 
nfaces=length (node3) ; 

%  load  vind  array  which  gives  the  three  nodes  for  each 
triangle 
for  i=l:nfaces 

pts= [nodel (i)  node2(i)  node3(i)]; 
vind (i, : ) =pts; 
end 

figured) 

%  This  section  plots  a  mesh 
if  icur==0 
for  i=l:nfaces 

X=[xpts (vind(i, 1) )  xpts (vind(i, 2) )  xpts (vind (i, 3)  ) 
xpts (vind (i, 1 ) ) ] ; 


Y=[ypts (vind(i, 1) )  ypts (vind (i, 2 ) 
ypts (vind(i,  1) )  ] ; 

Z=[zpts (vind(i, 1) )  zpts (vind(i,2) 
zpts (vind (i,  1) )  ] ; 

plots (X,Y,Z, 'm' ) 
if  i  ==  1 
axis  equal 
view (60,45) 
hold  on 
end 
end 
end 

hold  off 


)  ypts (vind (i,  3)  ) 
)  zpts (vind (i,  3 )  ) 
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%  FACT.M 

%  Subroutine  for  calculating  the  factorial  of  a  number 


function  f=f act (m) 

f=l; 

if  m>=l 

for  n=l:m 

f=f*n; 

end 

end 
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%  GF.M 

function  gf=gf(n,w) 

gf=  (exp  ( j*w)  -1)  /  (j*w)  ; 
if  n>=0 
for  in=l:n 

gf= (exp ( j  *w) -n*gf ) / ( j  *w) ; 
end 

end 
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%  RES.M 

%  Program  to  Calculate  Resistivity  on  Each  Triangle  using  a 
%  Linear  Taper 
for  n=l:nfaces 

xm  (n)  =  (r (max (vind (n, : ) ) , 1) +r (min (vind (n, : ) ) , 1) ) /2; 
ym (n)  =  (r (max (vind (n, : ) ) , 2) +r (min (vind (n,  : ) ) ,  2) ) /2; 
zm  (n)  =  (r (max (vind (n, : ) ) , 3) +r (min (vind (n, : ) ) , 3) ) /2; 

%  vector  to  midpoint  of  each  triangle 
mn (n, : ) = [xm (n)  ym(n)  zm(n)]; 

%  angle  from  z  axis  to  midpoint  vector 
if  zm(n)~=0 

omega  (n)  =pi/2-acos  (zm(n)  /  (sqrt  (xm(n)  ''2+ym(n)  '^2)  +le-8)  )  ; 
%arclength  from  z  axis  to  midpoint  of  each  triangle 
arc  (n)  =xm(n)  "'2+ym(n)  ''2; 
end 

%arclength  from  z  axis  to  midpoint  if  plate  is  not  flat 

arc (n) =omega (n) / (2*pi) ; 

end  %  end  of  triangle  loop 

rss=377/max (abs (arc) } ; 

end 
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%  SYNTH. M 

%  Resistivity  Synthesis-Monostatic 
clear  all 
format  long 

%  illumination  flag:  =0  external  face  only 
iflag=0; 

Lt=.001; 

Nt=2  ; 
wave=l ; 
bk=2*pi/wave; 
rad=pi/180; 

%  Incindent  wave  polarization 
Et=l+j*0;  %TM-z 

Ep=0+j*0;  %TE-z 

%  Load  results  from  potcos  program 
load  Ethr.m; 
load  Ethi.m; 

Ethii=Ethi*i; 

Eth=Ethr+Ethii ; 
load  Ephr.m; 
load  Ephi.m; 

Ephii=Ephi*i; 

Eph=Ephr+Ephii ; 
load  rs.m; 

%  Wave  ampltidue  at  all  vertices 

Co=l .  ; 

build 

ntria=nfaces; 
nvert=nvert s ; 
x=xpts; 
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y=ypts; 

z=zpts; 

title (' triangular  model  of  the  scattering  surface') 

xlabel ( 'x' ) 

ylabel ( ' y ' ) 

zlabel('z') 

for  i=l:nvert 

text (x  (i) -max (x) /20, y (i) -max (x)/20,z(i), num2str (i) ) 

end 

hold  off 

dispChit  return  to  continue  calculation') 
pause 
%  label 

%  Define  position  vectors  to  vertices 
for  i=l:nvert; 

r(i,:)  =  [x(i)  y(i)  z(i)]; 

end 

%  Get  edge  vectors  and  normals  from  edge  cross  products 
for  i=l:ntria 

A=r (vind(i,2) , : )-r (vind(i,  1)  ,  : )  ; 

B=r (vind (i, 3) , : ) -r (vind (i, 1) , ; 

C=r (vind(i, 3) , : ) -r (vind(i, 2) , : ) ; 

N (i, : ) =cross  (A, B) ; 

%  check  for  proper  direction  of  normal  vector 
if  N(i,3)<0 

N(i,  : ) =-N(i,  : )  ; 

end 

%  Edge  lengths  for  triangle  "i" 
d (i, 1) =norm (A) ; 
d (i, 2) =norm (B) ; 
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d (i, 3) =norm (C) ; 
ss=.5*suni(d(i,  : )  )  ; 

Area (i) =sqrt (ss* (ss-d (i,  l))*(ss-d(i,2))* (ss-d (i, 3) ) ) ; 
Nn=norm (N (i,  : ) )  ; 

N{i,  :)=N(i,  :)/Nn; 

Itransform  triangles  into  local  coordinates 

xl=[l  0  0] ; 

yl=[0  1  0]; 

zl=[0  0  1]  ; 

end 

%  Pattern  Loop 
it=0; 
k=  On¬ 
load  du.ra; 
for  ul=-l :du: 1-du 

for  vl=-l;du: 1-du; 

it=it+l ; 

%test  to  see  if  center  of  box  is  within  visible  region 
test- ( (ul+du/2) ^2+ (vl+du/2) ^2) ; 
if  test<=l 
k=k+l; 

phi=atan2 ( {vl+du/2)  ,  (ul+du/2) ) ; 
cp=cos (phi) ;  sp=sin (phi) ; 
theta=asin (sqrt (test) )  ; 
st=sin (theta) ;  ct=cos (theta) ; 

%  Center  of  box  is  used  so  u  and  v  are  redefined  using  center 
of  %  box 

u(k)=st*cp;  v(k)=st*sp;  w=ct; 
up=ct*cp;  vp=ct*sp;  wp=-st; 

%  Spherical  coordinate  system  unit  vectors 
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R(it, : ) = [u (k)  v(k)  w] ; 

%  Change  Efield  into  global  rectangular  coordinates 

strtra (1, : ) = [st*cp  ct*cp  -sp]; 

strtra (2, : ) = [st*sp  ct*sp  cp] ; 

strtra (3, : ) = [ct  -st  0] ; 

sei=  [0  Et  Ep] ; 

ei=strtra*sei ' ; 

%  Begin  loop  over  triangles 
siamt=0/ 
suiap=0; 
for  m=l:ntria 
t (m, : )=cross (zl,N (m, : ) ) ; 
beta (m) =acos (N (m, 3 ) ) ; 

alpha  (la)  =atan2  (t  (m,  2)  ,  t  (m,  1)  +le-6)  -pi/2; 

%  Test  to  see  if  front  face  is  illuminated 
ndotk=N (m, : ) *R (it, : ) ' ; 
if  ndotk  >=  0  &  iflag==0  |  iflag~=0 
%translation  into  local  coordinates 
%first  rotate  about  z  axis 

tral  (1, :)  =  [cos (alpha (m) )  sin (alpha (m) )  0]/ 
tral (2, : )  =  [-sin (alpha (m) )  cos (alpha (m)  )  0] ; 
tral (3, : )=[0  0  1] ; 

%  Rotation  about  yp  axis 

tra2  (1,  :)  =  [cos  (beta  (m)  )  0.  -sin  (beta  (m)  )]  ; 
tra2  (2,  :  )  =  [0  1  0]  ; 

tra2(3, ;)=[  sin (beta (m))  0  cos (beta (m) )] ; 
Dl=tral*R(it,  : ) ' / 

D2=tra2*Dl; 
upp=D2 ( 1 ) ; 
vpp=D2 (2) ; 
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wpp=D2 ( 3 ) ; 

%  Find  spherical  angles  in  local  coordinates 
stpp=sqrt  (upp''2+vpp''2 )  *sign  (wpp)  ; 
ctpp=sqrt  (l-stpp''2)  ; 
thpp=acos (ctpp) ; 
phipp=atan2 (vpp,upp+le-10) ; 

%  Phase  at  the  three  vertices  of  triangle  m;  monostatic  RCS 
needs  %  "2" 

Dp=2*bk* ( (x (vind (m, 1 ) ) -x (vind (m, 3) ) ) *u (k) + . . . 

(y (vind(m, 1) ) -y (vind(m, 3)))*v(k)+... 

( z (vind (m, 1 ) ) -z (vind (m,  3) ) ) *w) ; 

Dq=2*bk* ( (x (vind (m, 2 ) ) -x (vind (m, 3 ) ) ) *u (k) + . . . 

(y (vind (m, 2) ) -y (vind (m, 3)))*v(k)+... 

(z (vind (m, 2) ) -z (vind (m, 3) ) ) *w) ; 

Do=2*bk*  (x  (vind  (m,  3)  )  *u  (k)  +y  (vind  (m,  3)  )  *v  (k)  +z  (vind  (m,  3) )  *w)  ; 

%incident  field  in  local  rectangular  coordinates 

eip=tral*ei; 

eipp=tra2*eip; 

%  translation  of  E-field  into  local  spherical  coordinates 
strtra2= [stpp*cos (phipp)  ctpp*cos (phipp)  -sin(phipp); 

stpp*sin (phipp)  ctpp*sin (phipp)  cos (phipp) ;  ctpp  -stpp  0]; 
epp=strtra2''  (-1)  *eipp; 

%  Area  integral  for  general  case 
DD=Dq-Dp ; 
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expDo=exp ( j *Do)  ; 
expDp=exp ( j*Dp) ; 
expDq^exp ( j *Dq) ; 

%  Special  Case  1 
if  abs{Dp)  <  Lt  &  abs(Dq)>=  Lt 
sic=0. ; 
for  n=0:Nt 

sic=sic+ ( j*Dp) ^n/fact (n) * (-Co/ (n+1) +expDq* (Co*gf (n, -Dq) ) ) 
end 

Ic (k,m) =sic*2*Area (m) *expDo/j/Dq; 

%  Special  Case  2 

elseif  abs(Dp)  <  Lt  &  abs (Dq)  <  Lt 
sic=0 . ; 
for  n=0:Nt 

for  nn=0:Nt 

sic=sic+  ( j*Dp)  ^n*  ( j*Dq)  ''nn/fact  (nn+n+2)  *Co; 

end 

end 

Ic (k,m) =sic*2*Area (m) *expDo; 

%  Special  Case  3 

elseif  abs (Dp)  >=  Lt  &  abs (Dq)  <  Lt 
sic=0 . ; 
for  n=0:Nt 

sic=sic+  ( j*Dq)  ''n/fact  (n)  *Co*gf  (n+1,  -Dp)  /  (n+1)  ; 
end 

Ic (k,m) =sic*2*Area (m) *expDo*expDp; 

%  Special  Case  4 
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elseif  abs (Dp)  >=  Lt  &  abs (Dq)  >=  Lt  &  abs (DD)  <  Lt 
sic=0 . ; 
for  n=0:Nt 

sic  =  sic+ ( j*DD)  ^'n/fact  (n)  *  (-Co*gf  (n,  Dq) +expDq*Co/ (n+1)  )  ; 
end 

Ic (k,m) =sic*2*Area (m) *expDo/j/Dq; 
else 

Ic (k,m)  =2*Area (m) *expDo* (expDp*Co/Dp/DD-expDq*Co/Dq/DD-Co/Dp 
/Dq)  ; 

end 

%  end  of  special  cases  test 
end  %end  of  illumination  test 
end  %end  of  triangle  loop 
Eth2  (k)=Eth(it)  ; 

end  %  end  of  test  for  visible  region 

end  %  end  of  v  loop 

end  %  end  of  u  loop 

ac=pinv (Ic) *Eth2 ' ; 

for  k=l:ntria 

if  imag (ac (k) ) <=. 0001 
ac (k) =real (ac (k)  )  ; 
end 

sr (k)=(l/ac(k)-l) *.5; 

%  correct  round-off 
s (k) =round (sr (k) *le4 )  ; 
sr(k)=s (k)*le-4; 
end 

k=l :ntria; 
figure (2 ) 

plot (k, rs, k/ sr, ' - ' ) ; 
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title (' Comparison  of  Original  and  Synthesized  Resistivity') 
xlabel (' Triangle  Number') 
ylabel { 'Normalized  Resistivity') 
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%  POTESS.M 

%  PO  scattering  from  triangles 
clear  all 

%  illumination  flag:  =0  external  face  only 
iflag=l; 

Lt=.001; 

Nt=2  ; 
wave=l; 
bk=2*pi/wave; 
rad=pi/180; 

%  Incindent  wave  polarization 
Et=l+j*0;  %TM-z 

Ep=0+j*0;  %TE-z 

phi=0; 

phr=phi*rad; 

cp=cos (phr) ; sp=sin (phr) ; 

%  Wave  ampltidue  at  all  vertices 

Co=l . ; 

build 

ntria=nfaces; 
nvert=nverts ; 
x=xpts; 
y=ypts; 
z=zpts; 

title (' triangular  model  of  the  scattering  surface') 

xlabel { 'x' ) 

ylabel('y') 

zlabel('z') 

for  i=l invert 

text (x (i) -max (x) /20, y (i) -max (x) /20,  z (i) ,num2str  (i) ) 
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end 

hold  off 

disp('hit  return  to  continue  calculation') 
pause 

%  Define  position  vectors  to  vertices 
for  i=l;nvert; 

r (i,  : )  =  [x  (i)  y  (i)  z (i)  ]  ; 

end 

%Use  synthesized  resistivity 

%load  srO.m  %  0  resistivity 

%load  srl.m  %  377  ohm/sq  resistivity 

load  srl.m  %  linear  resistivity 

rs=srl; 

%  edge  vectors  and  normals  from  edge  cross  products 
for  i=l:ntria 

A=r (vind(i/ 2) , :)-r {vind(i, 1) ,  : ) ; 

B=r (vind (i, 3) , : ) -r (vind (i,  1) ,  : ) ; 

C=r (vind(i, 3) , : ) -r (vind(i,  2)  ,  : )  ; 

N{i, : ) =cross (A,  B) ; 

%  check  for  proper  direction  of  normal  vector 
if  N(i/ 3) <0 

N(i,  ;  )=-N(i,  : )  ; 

end 


%  Edge  lengths  for  triangle  "i" 
d(i,  1) =norm(A)  ; 
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d (i,  2)  =norin(B) ; 
d(i, 3) =norm(C) ; 
ss=.5*sum(d(i, : ) )  ; 

Area (i) =sqrt (ss* (ss-d (i, 1) ) * (ss-d (i, 2) ) * (ss-d (i, 3) ) ) ; 
Nn=norin (N  (i,  : )  ) ; 

N(i,  : )  =N(i,  : )  /Nn; 

%transform  triangles  into  local  coordinates 

xl=[l  0  0]; 

yl=[0  10]; 

zl=[0  0  1]; 

end 

%  Pattern  Loop 
start=0; 
stop=180; 
del=l; 

it=floor ( (stop-start) /del) +1; 
for  i=l:it 

theta (i) =start+ (i-1) *del; 
thr=theta (i) *rad; 
st=sin (thr) ;  ct=cos (thr)  ; 
u=st*cp;  v=st*sp;  w=ct; 
up=ct*cp/  vp=ct*sp;  wp=-st; 

%  Spherical  coordinate  system  unit  vectors 
R(i,  :  )  =  [u  V  w]  ; 

%  Change  Efield  into  global  rectangular  coordinates 

strtra (1, : ) = [st*cp  ct*cp  -sp] ; 

strtra (2,  : )  =  [st*sp  ct*sp  cp] ; 

strtra (3, :)= [ct  -st  0]; 

sei= [0  Et  Ep] / 
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ei=strtra*sei ' ; 

%  correct  orientation  for  normal  after  90cieg 
if  theta (i) ==91 
N=-N; 

end 

%  Begin  loop  over  triangles 
sumt=0 ; 
siimp=0 ; 
for  m=l:ntria 
t (m, : ) =cross (zl,N (m, : ) ) ; 
beta (m) =acos (N (m,  3 ) ) ; 

alpha (m) =atan2  (t (m, 2) , t (m, 1) +le-6) -pi/2; 

%  Test  to  see  if  front  face  is  illuminated 
ndotk=N (m, : ) *R(i, ; ) ' ; 
if  ndotk  >=  0  &  iflag==0  |  iflag~=0 
%translation  into  local  coordinates 
%first  rotate  about  z  axis 
tral (1/ ;)= [cos (alpha (m) )  sin (alpha (m) )  0]; 
tral (2, : ) = [-sin (alpha (m) )  cos (alpha (m) )  0]/ 

tral  (3,  :  )  =  [0  0  1]; 

%  Rotation  about  yp  axis 

tra2 (1, :)= [cos (beta (m) )  0  -sin (beta (m) )] ; 
tra2  (2,  :  )  =  [0  1  0]  ; 

tra2(3,:)=[  sin (beta (m))  0  cos (beta (m) )] ; 
Dl=tral*R (i, : ) ' ; 

D2=tra2*Dl; 
upp=D2 ( 1 ) ; 
vpp=D2 (2) ; 
wpp=D2 ( 3 )  ; 

%  Find  spherical  angles  in  local  coordinates 
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stpp=sqrt  (upp''2+vpp^2)  *sign  (wpp) ; 
ctpp=sqrt  ( l-stpp''2 )  ; 
thpp=acos (ctpp) ; 
phipp=atan2 (vpp, upp+le-10) ; 

%  Phase  at  the  three  vertices  of  triangle  m;  monostatic  RCS 
needs  %  "2" 

Dp=2*bk;*  (  (x  ( vind  (m,  1 )  )  -x  (vind  (m,  3 )  )  )  *u+  .  .  . 

(y (vind(m, 1) ) -y (vind(m, 3) ) ) *v+. . . 

(z  (vind(m, 1) ) -z (vind(m, 3) ) ) *w) ; 

Dq=2*bk* ( (x (vind (m, 2 ) ) -x (vind (m, 3 ) ) ) *u+ . . . 

(y (vind (m, 2) ) -y (vind (m, 3) ) ) *v+. . . 

(z (vind (m, 2) ) -z (vind(m, 3) ) ) *w) ; 

Do=2*bk* (x (vind (m,  3) ) *u+y (vind(m,  3) ) *v+z (vind(m, 3) ) *w) ; 
%incident  field  in  local  rectangular  coordinates 
eip=tral*ei; 
eipp=tra2*eip; 

%  translation  of  E-field  into  local  spherical  coordinates 
strtra2= [stpp*cos (phipp)  ctpp*cos (phipp)  -sin(phipp); 
stpp*sin (phipp)  ctpp* sin (phipp)  cos (phipp) ;  ctpp  -stpp  0] ; 
epp=strtra2''  (-1)  *eipp; 

%Reflection  coefficients 
gamp=0; 

if  (2*rs (m) +ctpp) ~=0 
gamp=l/ (2*rs (m) +ctpp) ; 
end 

gamn=l/ (2*rs (m) *cos (thpp) +1) ; 
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%Surface  Currents 

Jxpp  (i,in)  =  (epp  (2)  *gamp-epp  (3)  *gamn)  *ctpp^2; 

Jypp  (i,in)  =  (epp  (2)  *gamp+epp  (3)  *gainn)  *ctpp; 

%  Area  integral  for  general  case 

DD=Dq-Dp; 

expDo=exp(j*Do) ; 

expDp=exp ( j  *  Dp ) ; 

expDq=exp ( j*Dq) ; 

%  Special  Case  1 

if  abs (Dp)  <  Lt  &  abs(Dq)>=  Lt 

%disp('case  1') 

icase (i,m) -1; 

sic=0.; 

for  n=0:Nt 

sic=sic+  ( j*Dp)  ''n/fact  (n)  *  (-Co/  (n+1)  +expDq*  (Co*gf  (n,  -Dq)  ) ) ; 
end 

Ic  (i,]ii)  =sic*2*Area  (m)  *expDo/j  /Dq; 

%  Special  Case  2 

elseif  abs (Dp)  <  Lt  &  abs (Dq)  <  Lt 
%disp ( ' case  2 ' ) 

icase (i,m) =2; 
sic=0 . ; 
for  n=0:Nt 

for  nn=0:Nt 

sic=sic+  ( j*Dp)  ''n*  (j*Dq)  ^'nn/fact  (nn+n+2)  *Co; 
end 
end 

Ic  (i,in)  =sic*2*Area  (m)  *expDo; 
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%  Special  Case  3 

elseif  abs (Dp)  >=  Lt  &  abs (Dq)  <  Lt 
%disp('case  3') 

lease  (i,ia)  =3; 
sic=0 . ; 
for  n=0:Nt 

sic=sic+ ( j *Dq) ^n/fact (n) *Co*gf (n+1, -Dp) / (n+1 )  ; 
end 

Ic (i/m) =sic*2*Area (m) *expDo*expDp; 

%  Special  Case  4 

elseif  abs (Dp)  >=  Lt  &  abs (Dq)  >=  Lt  &  abs(DD)  <  Lt 
%disp('case  4') 

lease (i,m) =4; 
sic=0.; 
for  n=0:Nt 

sic  =  sic+ ( j*DD) ''n/fact  (n)  *  (-Co*gf  (n,  Dq) +expDq*Co/ (n+1)  ) ; 
end 

Ic (i,m) =sic*2*Area (m) *expDo/j/Dq; 
else 

lease (i,m) =0; 

Ic (i,m) =2*Area (m) *expDo* (expDp*Co/Dp/DD-expDq*Co/Dq/DD-Co/Dp 
/Dq)  ; 

end 

%  end  of  special  cases  test 

%Scattered  field  in  local  rectangular  coordinates 
Esxpp=Jxpp (i, m) *Ic (i, m) ; 

Esypp=Jypp (i,m) *Ic (i,m) ; 

Eszpp=0; 

%retranslation  to  global  rectangular  coordinates 
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Espp (m, : ) = [Esxpp  Esypp  Eszpp] ; 

Esl=tra2'"  (-1)  *Espp  (m,  : )  '  ; 

Es=tral^ (-1) *Esl; 

%retranslation  to  global  spherical  coordinates 
%seipp (m, : ) =strtra^ (-1) *Es; 
seipp (m, 1) =0; 

seipp (m, 2) =ct*cp*Es (1 ) +ct*sp*Es (2) -st*Es (3)  ; 
seipp (m, 3)=-sp*Es (l)+cp*Es (2) ; 

%sum  over  all  triangles  to  get  scattered  field 

sumt=sumt+seipp (m, 2 ) ; 

suinp=sump+seipp  (m,  3)  ; 

end  %end  of  illumination  test 

end  %end  of  triangle  loop 

Eth=sumt; 

Sth(i)=10*logl0  (4*pi*abs  (Eth)  ''2+l.e-lO)  ; 
Eph=sump; 

Sph (i) =10*logl0 (4*pi*abs (Eph) ^2+1 .e-10)  ; 
end 

Sth  ( : )  =max  (Sth  ( : )  ,  -50)  ; 

Sph { : ) =max ( Sph ( : )  ,  -50 ) ; 

%load  rcsO.m 
%load  rcsl.m 
load  rcsl.m 
%plot  res 
figure (2) 

plot (theta, rcsl, theta, Sth, ' — ' ) , grid 
title { ' solid: original  dashed: synthesized' ) 
xlabel ( 'monostatic  angle,  theta  (deg)') 
ylabel  ( '  sigma/wavl''2  (dB)  ' ) 
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